Total recoil: the maximum kick from nonspinning black-hole binary inspiral 
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When unequal-mass black holes merge, the final black hole receives a "kick" due to the asymmetric 
loss of linear momentum in the gravitational radiation emitted during the merger. The magnitude 
of this kick has important astrophysical consequences. Recent breakthroughs in numerical relativity 
allow us to perform the largest parameter study undertaken to date in numerical simulations of 
binary black hole inspirals. We study non-spinning black-hole binaries with mass ratios from q — 
M 1 /M 2 = 1 to q = 0.25 (n = q/(l + q) 2 from 0.25 to 0.16). We accurately calculate the velocity 
of the kick to within 6%, and the final spin of the black holes to within 2%. A maximum kick of 
175.2 ± 11 kms" 1 is achieved for n = 0.195 ± 0.005. 

PACS numbers: 04.25.Dm, 04.30.Db, 95.30.Sf, 98.80.Jk 



Introduction. — Anisotropic emission of gravitational 
waves from the coalescence of black-hole binaries carries 
away linear momentum and thus imparts a recoil on the 
merged hole. This recoil, often referred to as a "kick" 
or "rocket effect" , has important consequences for vari- 
ous astrophysical scenarios. The displacement or ejection 
of black holes as a result of a black-hole merger not only 
leads to a population of interstellar and intergalactic mas- 
sive black holes, but also has severe repercussions on the 
formation of supermassive black holes and the structure 
of the host galaxies [H, [3, S B 01 • The demography of 
massive black holes is also relevant for the expected num- 
ber of sources for the space-based gravitational- wave de- 
tector LISA. The gravitational recoil might also manifest 
itself directly in astrophysical observations, such as the 
discovery of bright quasi-stellar objects without a host- 
galaxy (see [liSMEl) an d the distorted morphology of 
x -shaped radio sources @, 0, [ll|. Accurate recoil esti- 
mates are important for all of these astrophysical models. 

The strongest contribution to the kick is made during 
the plunge and merger of the black holes. In this regime 
nonlinear general relativistic effects preclude reliable an- 
alytic treatment, and even the most recent sophisticated 
analytic estimates [l2, 1^, 14, 15, which give a maxi- 
mum kick varying from 50 to 500 km s _ 1 , carry uncertain- 
ties of 25 to 50%. Numerical studies using full numerical 
relativity are necessary to accurately determine the to- 
tal recoil. In contrast to analytic calculations, numerical 
simulations contain only one physical approximation: the 
initial data are not exactly equivalent to an astrophysical 
inspiral process. However, the resulting errors decrease 
as the black holes are placed further apart. If the initial 
separation is large enough that its effect on the estimate 
of the kick or final angular momentum is minimal, then 
the result can be said to be accurate and free from any 
physical approximation. 

Recent breakthroughs in numerical relativity 17|, Il8l . 
HiL [20! ] have made possible long-term stable numerical 
evolutions of binary black hole systems for several or- 
bits through merger and ringdown [2H I22L I23L 24, 25, 



26l . |27| . Drastic improvements in computational effi- 
ciency achieved by a new generation of accurate finite- 
difference mesh- refinement codes, e.g. 25|, l28(, pave the 



way for the large parameter studies necessary to ex- 
plore the parameter space of binary black hole inspiral. 
Here we present the first such study, comprising roughly 
three dozen unequal-mass initial-data sets, to compute 
the gravitational recoil and the spin of the final black 
hole. Implications for gravitational-wave data analysis 
will be explored elsewhere. 

Early numerical estimates of the recoil suffered from 
low numerical accuracy, or initial black-hole separations 
that were too small. (Compare, for example, the value 
of 240±140 kms" 1 for q = 0.5 in j2§| and the estimate 
33 kms -1 reported in [24j for q = 0.85.) The first ac- 
curate numerical result was recently given by Baker, et 
al. [30(, who found a value of 101±15 kms -1 for q = 0.67. 

In this paper we present results from numerical simu- 
lations of unequal-mass nonspinning black-hole binaries 
with mass ratios q = 1.0 to q = 0.253. We estimate that 
the total error in the kick velocities that we quote is less 
than 6% (see below). We are thus able to improve on the 
numerical and physical accuracy achieved in [30], and, 
more importantly, our numerical simulations for the first 
time cover a range of mass ratios large enough to accu- 
rately determine the maximum kick resulting from non- 
spinning binaries. We calculate the maximum kick to be 
175.2 ± 11 kms" 1 . 



Numerical methods. — Numerical 
performed with the BAM code [H [H| 
implemented the "moving puncture" 



ing black-hole binaries [18|, [19|, l31|. 



simulations were 
in which we have 
method of evolv- 
This approach is 

based on initial data of puncture type 32, 33|, which are 
evolved using the BSSN/(3+l) formulation of the Ein- 
stein equations [3, l35j |- with a suitable gauge choice. 
The code uses a box-based mesh-refinement grid struc- 
ture with coarser refinement levels being centered on 
the origin and the high resolution levels consisting of 
two components centered around each hole and follow- 
ing their motion across the computational domain. The 
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evolution uses a fourth-order accurate Runge-Kutta inte- 
grator, and Berger-Oliger timestepping for the mesh re- 
finement. Gravitational waves are extracted in the form 
of the Newman-Penrose scalar ^4 on spheres of constant 
coordinate radius. Details of all aspects of the implemen- 
tation have been presented in [28j . 

In order to determine the physical parameters of the 
initial data, we estimate the initial momenta of the black 
holes using the 3PN-accurate formula given in Sec. VII 
in [28(, for given masses and coordinate separation. The 
black holes have masses Mi and M2, and the total black- 
hole mass is M — Mi + M%. The total gravitational 
energy of the system is Madm [36]. In the notation of 
[2^, we use the Xv=2 method, and label our runs by the 
number of gridpoints, i, on the finest level. We have 
observed that grid sizes of at least i = 48 are required 
to achieve fourth-order convergence in the waveforms and 
puncture trajectories. Here we have performed runs with 
i = 32,40,48,56 for all mass ratios, corresponding to 
finest resolutions of 1/26, 1/32, 1/38 and 1/45 (for de- 
tails see Table I in [28j). The main purpose of the low 
resolution runs was to develop and test our strategy for 
setting up our simulations. In order to confirm that our 
i = 56 runs are in the fourth-order convergent regime 
for the waveforms and puncture tracks, we have per- 
formed convergence tests using resolutions characterized 
by i = 56, 64, 72, 80 (with resolutions at the punctures of 
1/45, 1/51, 1/58, and 1/64) for selected mass ratios. We 
start with a simulation of an equal-mass binary, that is, 
each black hole has an initial mass of Mi = M2 =0.5. 
The mass ratio is then changed by increasing the mass 
of one of the holes. For example, the black holes have 
masses Mi — 0.5 and M2 = 1.0 in the case of a mass- 
ratio q = M1/M2 = 0.5. The idea behind this strategy 
is to keep constant the effective numerical resolution of 
the small black hole, while increasing that for the larger 
hole. For consistency, it is then necessary to proportion- 
ately increase the initial distance between the black holes, 
the distance of the outer boundary of the numerical grid, 
and the extraction radii of the gravitational waves. We 
find that sufficiently accurate (for the purpose of this 
study) gravitational-wave signals could be extracted at 
30M. Finally, the radiated linear momentum is calcu- 
lated from ^4 according to 



dt 



lim 



16tt 



dn 



(1) 



37|. 



where li = (sin 6 cos <f>, sin 9 sin (/>, cos 9) [3C 

The puncture initial data with Bowen-York extrinsic 
curvature used in our simulations arc known to con- 
tain spurious gravitational radiation. This pulse quickly 
leaves the system, but carries away linear momentum 
which wc find to be in the direction transverse to the 
black-hole motion. This small initial "kick" of the order 
of lOkrns -1 is not part of the astrophysical black- hole 



recoil we wish to measure. We therefore wait for the ini- 
tial pulse to pass through the extraction sphere (after 
about 50 Mad m in our simulations) , and calculate Pi by 
starting the integration of dPi / dt as given in Eq. HJ at 
to = 50Madm ■ In the course of the inspiral, dPi / dt os- 
cillates around zero, and, thus, the final integrated kick 
Pi will depend on the choice of to. For the separations 
we consider ambiguities in this choice introduce an un- 
certainty in Pi of 3%. 

Results. — We performed two sets of runs. In the first 
set, the mass ratio was varied from q = 1.0 to q = 0.25, 
and the initial coordinate separation of the black holes 
was kept fixed at ro = 7.0M. Convergence tests were 
performed for mass ratios q — 0.4, 0.33, 0.286 with finest- 
grid resolutions of hi — 1/45, h% — 1/51 and I13 = 1/58. 
In order to test the convergence properties of the recoil, 
we study the components P x and P y of the final kick. 
We find the results to be consistent with second-order 
convergence, as shown for q = 0.33 in Fig. Q] The kick 
is calculated by integrating twice in time a function that 
is several orders of magnitude smaller than the wave sig- 
nal ^4(2, 2), and we believe this to be the reason for the 
lower-order accuracy than observed in [28| . Higher res- 
olutions, currently beyond our computational resources, 
are necessary to achieve the limiting overall convergence 
behavior of the code for all quantities. There is a large 
error in the estimate of the kick's direction, but the kick's 
magnitude is calculated with high accuracy. We estimate 
that the total kicks are calculated with a numerical error 
of less than 2%. 

From Fig. Q] we see that v x — P x /Mp and v y — P y /Mp, 
where Mp is the mass of the final black hole, oscillate 
around zero during the inspiral. This is consistent with 
a small continuous loss of linear momentum, like water 
from a spinning lawn sprinkler, which pushes the center- 
of-mass in a roughly circular motion, before the final kick 
during merger. 

For the second set of runs we fix rj = 0.19, but vary the 
initial separations, setting r — 6. OA/, 7.0M and 8.0M. 
These simulations demonstrate the extra contribution to 
the kick at larger ro to be small, so that we estimate the 
accumulated kick from the early inspiral to be less than 
1%. This is comparable to post-Newtonian estimates of 
the contribution to the kick up to our initial separation 

Figure[2]shows the total kick from all of these runs, as a 
function of 77. Combining all the errors we have discussed 
(3% due to the choice of to in Eqn. (JJ), 2% numerical 
error, and 1% due to the neglected contribution from 
earlier inspiral), we conservatively estimate a total error 
of less than 6%. A least-squares fit of the kick with v = 
Ar] 2 ^l — 477(1 + Brf) (based on the formula of Fitchett 
Q) gives ,4= 1.20xl0 4 and B = -0.93 with X 2 = 48.25 
with v — 29 degrees of freedom. From our curve fit we 
calculate a maximum kick of V max = 175.2 ± llkms -1 
at 77 = 0.195 ± 0.005 (q = 0.36 ± 0.03). This agrees 
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FIG. 1: Components v x and v y of the kick velocity as functior 
of time, for r/ = 0.19, for resolutions hi = 1/45, hi = 1/51 
and /13 = 1/58. Top panels: kick components v x and v y . 
Lower panels: demonstration of second-order convergence. 



FIG. 3: Total kick velocity (v — y/vE+Vy) as function of 
time for r\ = 0.19. Top panel: for resolutions hi, hi and hz as 
described in the text. Lower panel: for runs with three initial 
black-hole separations r = G.0M, 7.0M, 8.0M. 



with the estimate of 114 ± 65 kms -1 of 14|, and also 
the close-limit analyses in l]| 16 1. The higher estimate 
of 250 ± 50km/s reported in Ref. [l]| does not include 
a possible "breaking" effect in the ringdown phase, and 
their values agree well with the local maximum in Fig. [3] 
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FIG. 2: The kick velocity as a function of mass ratio, with an 
error of ±6% indicated by the dotted lines. We also indicate 
previous numerical results from Baker, et al [30j| . Campanelli 
[291 |. and Herrmann, et a/ 1241 . and the analytic estimates of 
Damour and Gopakumar [l4| and Sopuerta, et al [TBI ]. 



Finally, we address the spin of the merged black hole. 
An understanding of the demographics of black holes, 
in particular the expected values of spins, is of essential 
importance for astrophysics, and also for developing ap- 
proaches to explore the astrophysically relevant binary 




FIG. 4: Left panel: Radiated energy and angular momentum 
as function of the mass ratio. Right panel: Spin of the final 
black hole (a/A//) as function of the mass ratio. The spin 
curves can be fit by a/M f = 0.089(±0.003) + 2.4(±0.025)??. 



black hole inspiral parameter space. We have computed 
the initial angular momentum from surface integrals at 
the wave extraction sphere as described in [5§|, and the 
final angular momentum from computing the wave ring- 
down frequency from an amplitude-phase decomposition 
of the radiation signal, and comparing with the depen- 
dence of angular momentum of a Kerr black hole on the 
ringdown frequency as quoted in [39] . An error estimate 
is obtained from evaluating the angular momentum sur- 
face integrals at the extraction sphere at the end of the 
simulation. We thus find our results to be accurate to 
within about 2 %. Our results for the unequal-mass se- 
quence considered here are displayed in Fig. [4] In the 
regime we consider, the dependence of the final spin on 
the mass ratio is approximately linear when expressed as 
a function of rj: a/M f = 0.089(±0.003) + 2.4(±0.025)7? 
(the correct result for rj = corresponds to a = 0). 
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Discussion. — Our results constitute the first accu- 
rate numerical calculations of the recoil velocity from the 
merger of nonspinning black holes for a large range of 
mass ratios, q = 1.0 to 0.253. This range allows us to de- 
termine the maximum kick as V max — 175.2 ± 11 kms -1 
at a mass ratio of q = 0.36 ± 0.03. Theoretical esti- 
mates predict a maximum kick at q — 0.38, and agree 
well with our value. Our maximum kick velocity is within 
the error bounds of the analytic estimates of Damour and 
Gopakumar [3] and Sopuerta, et al fl5| . which have sig- 
nificantly refined an earlier estimate by Favata et al 12 1 . 
The kick before ringdown agrees well with the results in 
[l3j |. We also find consistency with the numerically cal- 
culated value 105 ± 10 kms" 1 for q = 0.67 in 

The results presented here required the largest param- 
eter study of numerical binary black hole evolutions to 
date, comprising approximately three dozen data sets 
describing unequal-mass nonspinning inspiraling black 
holes. The efficiency of our code |28[ has enabled us to 
perform all the simulations quoted here at a total com- 
putational cost (including the development of our strat- 
egy and setup, production runs and convergence tests) of 
about 150 000 CPU hours. 

In future work we plan to study the consequences of our 
results for gravitational-wave data analysis. Extending 
this work to much higher mass ratios would require larger 
evolution times and higher resolutions, making such sim- 
ulations computationally very expensive, but in principle 
possible with our current techniques. The situation is 
similar if we wish to consider larger initial separations. 

We thank Emanuele Berti for giving us access to data 
for the quasinormal mode frequencies of Kerr spacetime, 
Clifford Will and Luc Blanchet for discussions of the ring- 
down breaking, and Iris Christadler at the Leibniz Com- 
puter Center (LRZ) for assistance in code optimization. 
Computations where performed at HLRS (Stuttgart) and 
LRZ (Munich). This work was supported in part by 
DFG grant SFB/Transregio 7 "Gravitational Wave As- 
tronomy" . 
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